# =====================================================================
# Supplementary File 5
# Analysis code for: Intra-wound vancomycin powder for the prophylaxis
# of surgical site infections in spinal surgery: a systematic review
# and meta-analysis
#
# Reproduces every pooled estimate, heterogeneity statistic, sensitivity
# analysis and publication-bias test reported in the manuscript.
#
# Requires: meta (>= 6.0), metafor (>= 4.0)
#   install.packages(c("meta", "metafor"))
# Analytical data: reproduced from Supplementary File 4 within this script
# =====================================================================

library(meta)
library(metafor)

## ---------------------------------------------------------------- data
## event.e / n.e = vancomycin powder arm; event.c / n.c = control arm

overall <- data.frame(
  study   = c("O'Neill 2011","Tubaki 2013","Hill 2014","Chotai 2017","Hey 2017",
              "Devin 2018","Mirzashahi 2018","Oktay 2021","Ushirozako 2021",
              "Vakayil 2021","Salimi 2022","Khalid 2023","Sonbol 2023"),
  event.e = c(0, 7, 5, 20, 1, 21, 10, 2, 21, 10, 13, 1, 1),
  n.e     = c(56, 433, 150, 1215, 117, 966, 193, 102, 623, 221, 187, 39, 126),
  event.c = c(7, 8, 11, 40, 17, 56, 5, 7, 33, 5, 25, 5, 6),
  n.c     = c(54, 474, 150, 1587, 272, 1090, 187, 107, 638, 221, 188, 39, 209),
  design  = c("NRS","RCT","NRS","NRS","NRS","NRS","RCT","NRS","NRS",
              "NRS","RCT","NRS","RCT"),
  cdc     = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE,
              TRUE, FALSE, FALSE, TRUE),
  stringsAsFactors = FALSE)

deep <- data.frame(
  study   = c("O'Neill 2011","Tubaki 2013","Hill 2014","Chotai 2017","Hey 2017",
              "Devin 2018","Adhikari 2020","Oktay 2021","Ushirozako 2021",
              "Vakayil 2021","Salimi 2022"),
  event.e = c(0, 6, 1, 20, 1, 7, 2, 1, 12, 3, 7),
  n.e     = c(56, 433, 150, 1215, 117, 966, 88, 102, 623, 221, 187),
  event.c = c(5, 6, 6, 40, 10, 42, 1, 4, 18, 1, 6),
  n.c     = c(54, 474, 150, 1587, 272, 1090, 70, 107, 638, 221, 188),
  stringsAsFactors = FALSE)

superficial <- data.frame(
  study   = c("O'Neill 2011","Tubaki 2013","Hill 2014","Hey 2017","Oktay 2021",
              "Ushirozako 2021","Salimi 2022"),
  event.e = c(0, 1, 4, 0, 1, 9, 5),
  n.e     = c(56, 433, 150, 117, 102, 623, 187),
  event.c = c(2, 2, 5, 7, 3, 15, 7),
  n.c     = c(54, 474, 150, 272, 107, 638, 188),
  stringsAsFactors = FALSE)

## ------------------------------------------------------- helper function
## Mantel-Haenszel fixed effect is the pre-specified primary model.
## DerSimonian-Laird random effects is reported alongside for every outcome.

run_ma <- function(d, label) {
  m <- metabin(event.e = d$event.e, n.e = d$n.e,
               event.c = d$event.c, n.c = d$n.c,
               studlab = d$study,
               sm = "OR", method = "MH", MH.exact = FALSE,
               method.tau = "DL", incr = 0.5, allincr = FALSE,
               common = TRUE, random = TRUE)
  cat("\n==== ", label, " (k = ", m$k, ") ====\n", sep = "")
  print(summary(m))
  invisible(m)
}

## -------------------------------------------------- primary and secondary
m_overall     <- run_ma(overall,     "Overall surgical site infection")
m_deep        <- run_ma(deep,        "Deep surgical site infection")
m_superficial <- run_ma(superficial, "Superficial surgical site infection")

## forest plots (Figures 3-5)
forest(m_overall,     common = TRUE, random = FALSE, xlim = c(0.005, 200))
forest(m_deep,        common = TRUE, random = FALSE, xlim = c(0.005, 200))
forest(m_superficial, common = TRUE, random = FALSE, xlim = c(0.005, 200))

## ---------------------------------------------------- sensitivity analyses
run_ma(subset(overall, design == "RCT"), "S1 Randomized controlled trials only")
run_ma(subset(overall, design == "NRS"), "S2 Non-randomized studies only")
run_ma(subset(overall, event.e > 0 & event.c > 0), "S3 Excluding zero-event studies")
run_ma(subset(overall, study != "Hey 2017"),
       "S4 Excluding the study at serious risk of bias")
run_ma(subset(overall, !(study %in% c("Hill 2014", "Chotai 2017"))),
       "S5 Excluding studies with baseline imbalance or unconfirmed univariable findings")
run_ma(subset(overall, cdc), "S6 Excluding studies without explicit CDC/NHSN definitions")
run_ma(subset(overall, study != "Chotai 2017"), "S7 Excluding the largest study")

## leave-one-out
cat("\n==== Leave-one-out, overall surgical site infection ====\n")
print(metainf(m_overall, pooled = "common"))

## ------------------------------------------------------- publication bias
cat("\n==== Funnel plot and Egger test (Figure 7) ====\n")
funnel(m_overall, common = TRUE, random = FALSE)
print(metabias(m_overall, method.bias = "linreg", k.min = 10))

sessionInfo()